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The two-dimensional Hubbard model defined for topological band structures exhibiting a quan- 
tum spin Hall effect poses fundamental challenges in terms of phenomenological characterization 
and microscopic classification. In the limit of infinite coupling U at half filling, the spin model 
Hamiltonians resulting from a strong coupling expansion show various forms of magnetic ordering 
phenomena depending on the underlying spin-orbit coupling terms. We investigate the infinite U 
limit of the Kane-Mele Hubbard model with z-axis intrinsic spin-orbit coupling as well as its gener- 
alization to a generically multi-directional spin orbit term which has been claimed to account for the 
physical scenario in monolayer Na2ir03. We find that the axial spin symmetry which is kept in the 
former but broken in the latter has a fundamental impact on the magnetic phase diagram as we vary 
the spin orbit coupling strength. While the Kane-Mele spin model shows a continuous evolution 
from conventional honeycomb Neel to XY antiferromagnetism which avoids the frustration imposed 
by the increased spin-orbit coupling, the multi-directional spin-orbit term induces a commensurate 
to incommensurate transition at intermediate coupling strength, and yields a complex spiral state 
with a 72 site unit cell in the limit of infinite spin-orbit coupling. From our findings, we conjecture 
that in the case of broken axial spin symmetry there is a large propensity for an additional phase 
at sufficiently large spin-orbit coupling and intermediate U. 

PACS numbers: 31.15.V-, 75.10.Jm, 03.65.Vf 



I. INTRODUCTION 



The discovery of the quantum Hall effect has initial- 
ized the era of topological phases in condensed matter 
physics. For non- interacting band structures with topo- 
logically unconventional properties, topological indices 
take over the role of conventional order parameters and 
can be linked to quantization phenomena of edge modes 
measured in experiment. The first example of such an in- 
dex was introduced by Thouless, Kohmoto, Nightingale, 
and den Nijs (TKNN) for the integer quantum Hall effect 
(IQHE)pJThey could show that the first Chern number — 
the TKNN invariant — is proportional to the transversal 
Hall conductivity a xy which is the integral of the Berry 
curvature over the Brillouin zone. Nearly a decade ago 
after Haldane realized that one can define lattice versions 
of IQHE called Chern insulators where complex hopping 
breaks time-reversal symmetry 2 } the most recent exam- 
ple of a non-interacting topological state of matter is the 
topological insulatop^. It is characterized by a Z2 topo- 
logical index^l Z2 topological insulators (TIs) have not 
only been proposed theoreticall}^^ but have also been 
found in subsequent experiments P The minimal model 
of a Z2 topological insulator is a four-band model pos- 
sessing a finite Z2 invariant, which in its simplest form 
is a minimal time-reversal invariant generalization of a 
Chern insulator. All two-dimensional band structures 
exhibiting a non-trivial Z2 invariant can be adiabatically 
transformed into each other, i.e. without closing the bulk 
gap. In contrast, transforming a Z2 TI phase into any 
other topologically trivial phase causes a quantum phase 
transition where the bulk gap must close. To date, these 



topological band insulators are well understood and sys- 
tematically classified by symmetry! 10 * 1 ^ 

As soon as interactions are taken into account, the 
full scope of possible scenarios extends to (i) topologi- 
cal band structure phases where the interactions would 
only renormalize the band parameters but do not change 
the topology along with (ii) conventional ordering phe- 
nomena where all features of the topologically non-trivial 
phase are go ne^^, and (iii) topological Mott insula- 
tors 5 E2Hl4125j anc j topological bulk order driven by 
strong interactions along with finite quant um di men- 
siorP^, fractionalization of quantum numberd 27 * 28 *, and 
fractional statistics 29 as well as interaction-driven topo- 
logical band structure phases which a re not a diabatically 
connected to the non-interacting limitP^^ 2 ^. In analogy 
to the non- interacting counterpart, topological bulk order 
was firs t disc overed in the fractional quantum Hall effect 
(FQHE) 28 30 before the concept of topological order was 
established by WerP^l Due to the diversity of possible 
phases even in the same symmetry sector, a general clas- 
sification for interacting topological phases is lacking so 
far. Aside from many other challenges, it is of particular 
interest whether the concept of a topological band struc- 
ture and topological bulk order can both manifest itself 
in a single microscopic modePH For example, competing 
magnetic fluctuations originating from a topological band 
structure model could manage to stabilize a topological 
spin liquid phase. This is the general motif of a class of 
scenarios which we further investigate in this article. 

As we will show in detail in the following, the pres- 
ence or absence of axial spin symmetry stemming from 
the topological band structure in the interacting case 





FIG. 1. Intrinsic spin orbit terms with amplitude A according 
(a) to Q for the Kane-Mele model and (b) to (2| for the SI 
model with multi-directional SOC amplitude A. 



will crucially determine the magnetic order and disorder 
phenomena which appear in the strong coupling limit. 
Generically, the full SU(2) is broken for interacting topo- 
logical band structure models because of spin orbit cou- 
pling terms. Still, it is both possible that the spin orbit 
terms break SU(2) down to U(l), leaving a continuous 
axial spin symmetry intact, or completely break spin ro- 
tation symmetry. Since its custodial time-reversal sym- 
metry is unaffected, it is irrelevant for the Z2 index of 
the weakly coupled model whether the axial spin sym- 
metry of the TI is conserved or not: although it has 
been shown recently that breaking of axial spin symme- 
try causes a momentum-dependent rotation of the spin 
quantization axis of the helical edge states p^J the topo- 
logical band structure with conserved spin symmetry can 
still be transformed into one with broken spin symmetry 
without closing of the bulk gap. In contrast, for strong 
interactions, the resulting phase diagram crucially de- 
pends on presence or absence of axial spin symmetry; 
more specifically, it was claimed that the combination of 
strong interactions and strong spin orbit coupling might 
give rise to a topologically ordered phase on the hon- 
eycomb lattice when spin is not conserved. This would 
then be a paradigmatic candidate model which includes 
both a topological band structure phase and topological 
bulk order in its phase diagram!^ Unfortunately, only 
the conserved U(l) symmetry appears to open up the 
possibility to successfully perform quantum Monte Carlo 
(QMC) simulations for the regime of intermediately cou- 
pled topological band structure models; when this sym- 
metry is absent, we instead have to rely on limited mean- 
field, slave-particle, or other approximate methods. 

In our work, we propose the strategy to first gain 
insight about this kind of models in the limit of in- 
finitely large interactions on the footing of an accurate 
method adapted to this limit, and to find out which 
of the approximate results at intermediate interaction 
strength is compatible with it. For this purpose, we 
employ pseudofermion functional renormalization group 
(PFFRG) which has been recently developed and em- 
ployed by two of us in the context of various models of 



frustrated magnet isnP^HMl, i n particular, the anisotropic 
spin terms do not pose additional challenges to the per- 
formance of the PFFRG, which at the same time allows 
us to study large system sizes beyond any other mi- 
croscopic numerical procedure for two-dimensional spin 
models. 

In this paper, we investigate the strong coupling limit 
of two different topological band structures accompanied 
with Hubbard onsite interactions on the honeycomb lat- 
tice: the Kane-Mele (KM) mo deP^ preserving axial spin 
symmetry and a related model which was proposed in the 
context of Na2lr03 by Shitade et a/P^ which explicitly 
breaks axial spin symmetry. Because of its connection to 
sodium iridate, it will be referred to as SI model in the 
following. We find that while magnetism in the presence 
of axial spin symmetry can generically avoid the frustra- 
tion effects caused by the anisotropic spin terms induced 
by spin-orbit coupling and generically yields commensu- 
rate magnetism, the broken axial spin symmetry scenario 
naturally leads to commensurate-incommensurate transi- 
tions and, as a consequence, a much more complex mag- 
netic phase diagram. As such, we conjecture that the 
latter scenario will be most promising to stabilize un- 
conventional, possibly topologically bulk ordered phases 
resulting from anisotropic spin terms. We also discuss 
our findings in the context of recent results^ for the cor- 
responding Hubbard models at finite coupling. 

The paper is organized as follows. In Section [TTJ we 
introduce the KM and SI models and discuss their main 
properties. The mean field phase diagrams of the cor- 
responding Hubbard models - the Kane-Mele-Hubbard 
(KMH) model as well as the sodium iridate Hubbard 
(SIH) model - are briefly reviewed in Section III We 
subsequently introduce the corresponding spin models in 
Section [IV| In Section |V| we elaborate on the PFFRG 
method which we employ to investigate the magnetic 
phase diagrams of the KM and SI spin models the results 



of which are presented in Section |VI| In Section |VII[ we 
draw a line from our findings at infinite coupling to the 
corresponding Hubbard models at finite coupling in the 
context of the recently proposed QSH* phase, a topolog- 
ically ordered phase in the SIH modelP^ In particular, 
we also point out important generalizations of our study 
with respect to Rashba coupling, which will generically 
break axial spin symmetry. In Section |VHI[ we conclude 
that the role of the axial spin symmetry is crucial to 
characterize magnetic order and disorder phenomena of 
interacting topological honeycomb band structures and 
leads to a better understanding of the general theme of 
interaction effects in topological insulators. 

Throughout this paper we use the following notations: 
the non-interacting topological insulators, i.e. the band 
structures are denoted by Hkm and h$i, respectively. The 
corresponding Hubbard models are called Hkm and Hsi 
while the spin models are denoted by Hkm and Hsi, re- 
spectively. The real nearest neighbor hopping amplitude 
is t; the intrinsic spin orbit couplings are called A for the 
KM model and A for the SI model. 
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II. TOPOLOGICAL BAND STRUCTURES 



B. Sodium iridate tight binding model 



The QSH honeycomb models are particularly accessi- 
ble from a theoretical perspective: as there are already 
two sites per unit cell, it is sufficient to study a single 
orbital scenario where complex hoppings generate the 
band inversion giving rise to a non-trivial Z2 invariant. 
There is hope that the QSH effect on the honeycomb lat- 
tice might be realized, e.g. by doping heavy adatoms in 
graphene 4 ^ or by using silicene 4 ^ which has recently been 
accomplished experimentally 4 ^. Depending on the con- 
cise form of the spin-orbit coupling terms, the axial spin 
symmetry may or may not be broken in the interacting 
case. In this section we briefly introduce the two repre- 
sentative models for both scenarios which are subject to 
further investigation in the following. 



A. Kane— Mele model 

Kane and Mele^ proposed the quantum spin Hall 
(QSH) effect in graphene based on symmetry considera- 
tion. They realized that a mass term cx a z T z rf does not 
violate any symmetries of graphene and thus must be al- 
lowed. Here, a is associated with the electron spin, r with 
the valleys, and 77 with the sublattices. The Kane-Mele 
model is governed by the tight-binding Hamiltonian 

^KM = ~t ^ 4a C ja + iX U ij C ia<(3 C j (3 (1) 

(ij)cr <ij> a(3 

In principle, there is also the Semenoff mass term which 
we will ignore for the moment. Similarly, the Rashba spin 
orbit term with amplitude Ar is neglected unless noted 
otherwise. The first term in ([!]) is the usual nearest- 
neighbor hopping on the honeycomb lattice giving rise 
to the Dirac band structure. The second term in is 
the lattice version of the <j z r z r] z -term (a second neigh- 
bor hopping) which corresponds to an intrinsic spin orbit 
coupling (SOC). The convention of this hopping is illus- 
trated in Fig.fT^i. The nearest neighbor hopping term 
preserves the Cq v lattice symmetry of the honeycomb 
lattice as well as SU(2) symmetry of the electron spin. 
The intrinsic SOC reduces the lattice symmetry to Cs v 
and the spin symmetry to U(l). Any finite A opens the 
gap of the Dirac band structure and gives rise to QSH 
effect, i.e. to a topological insulator phase characterized 
by a finite Z2 invariant, or, in this case, Chern number for 
each spin species. This situation is very special since the 
Hamiltonian fully decouples into two independent Chern 
insulators with opposite Hall conductivity. Generically, 
we expect the presence of additional terms breaking the 
U(l) spin symmetry and mixing the spin channels. The 
Rashba term is such an additional term which will be fur- 
ther commented on in Section lVIII Even for finite Rashba 
coupling A^, however, the QSH phase is stable as long as 
A i? <2 V / 3AP 



Soon after Kane and Mele's milestone works, it turned 
out that the spin orbit gap in graphene is vanishingly 
small. Therefore other materials with effective honey- 
comb structure were considered as candidates for the 
QSH effect as proposed by Kane and Mele. In 2008, Shi- 
tade et alW* came up with the sodium iridate Na2lr03 
as a layered honeycomb system. The authors claimed 
that the QSH effect might be realized if Coulomb inter- 
actions are not too strong. A monolayer was shown to 
be described by a Kane-Mele-type Hamiltonian. The in- 
trinsic spin orbit coupling was assumed to be relatively 
large due to the heavier iridium atoms in contrast to 
graphene's carbon atoms. Assuming trivial hybridiza- 
tion between nearest neighbor Ir atoms, Shitade et at. 
found an intrinsic SOC being similar but different to the 
KM SOC. It depends on the direction of the spin orbit 
hopping whether the spin degree of freedom is associated 
with <j x , (7^ , or a z . The sodium iridate model is governed 
by the Hamiltonian 

ft si = ~t c i* c j* + iX J2 H c \A c jp , ( 2 ) 

{ij)cr <ij> 7 a/3 

where 7 = x,y,z is associated with the different next- 
nearest neighbor links on the honeycomb lattice (Fig.J^)). 
The main difference of this generalized SOC compared 
to the KM SOC is that axial spin symmetry is not con- 
served. As for the KM model, infinitesimally small A 
opens the gap at the Dirac cones and causes QSH effect. 

The band structures of Hkm and hsi both belong to 
the Z2 universality class and are thus adiabatically con- 
nected. Both systems exhibit helical edge states on open 
geometries such as cylinder or disk. 

III. CORRELATED TOPOLOGICAL 
INSULATORS 

Let us now add Hubbard onsite interactions, 

Hi = uJ2 n it n a (3) 

i 

which yields rich phase diagrams for both band struc- 
tures. While the U-X phase di agram of the KMH model 
is well understood 113 1 15 1 16 1 19 1 43 1 44 1 , the U -X phase diagram 
of the SIH model is rarely studied 33 39 , and the available 
results are controversial. In the following, we will briefly 
review the phase diagrams of both Hubbard-type models. 

A. Kane— Mele— Hubbard model 

The KMH model is described by a combination of the 
KM and Hubbard model, 

Hkm = ^km + Hi . (4) 
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FIG. 2. (color online), (a) phase diagram of the Kane-Mele-Hubbard model as obtained in Ref. 1131 The transition from a 
topological insulator (Tl) to XY-plane antiferromagnet (AFM) was derived within slave-rotor theory which underestimates 
U c - (b) mean field phase diagram of the sodium-iridate Hubbard model as obtained in Ref.1331 The transition from Tl to a 
valence bond solid (VBS) phase that links to AFM was derived within slave-spin theory which overestimates U c - The phase 
diagram of (b) is qualitatively similar to (a) apart from the additional "QSH* phase". At A = A = and not too large U the 
semi-metal (SM) phase of graphene is present. See main text for details. 



In Ref.[13] the phase diagram shown in Fig.[2^i was de- 
rived through slave rotor theory. The semi-metal (SM) 
phase of graphene (A = 0) as well as the topological in- 
sulator phase (A ^ 0) are stable up to moderate inter- 
actions. Above a critical interaction strength U c , one 
finds an antiferromagnetically ordered phase which is of 
Neel type (A = 0) or of XY-type (A / 0), respectively. 
At A = and intermediate U, a quantum spin liquid 
phase has been proposed^ recently; this conjecture has 
been challenged lately. 46 For ve ry small A it survives but 
eventually vanishes for A < 0.05 g 15 l 16 l 43 l 44 l Since the spin 
liquid is destroyed by finite A and just a remnant of the 
non-topological A = case, we omit the phase here for 
clarity. Also for the strong coupling analysis in this paper 
we will assume that we are deep in the strong coupling 
regime where this intermediate coupling phenomenon is 
irrelevant for our analysis. 

B. Sodium iridate Hubbard model 

Recently, Ruegg and Fiete have studied the SIH 
model^ governed by the Hamiltonian 

Hsi = h S i + Hi • (5) 

They used a Z 2 slave-spin mean-field approach and pro- 
posed an interesting phase diagram (Fig.J^j)). It is simi- 
lar to the KMH model, while there is an additional phase 
for large SOC A and large J7, dubbed QSH* phase, which 
presumably extends to the strong coupling regime. Note 
that this is not a quantum spin Hall phase, but a topo- 
logical liquid which is characterized by a four-fold de- 
generacy on a torus, where the elementary excitations 
are fractional particles obeying Abelian statistics. Re- 
cently it was questioned, however, whether the employed 
Z2 slave spin approach is justified. 47 Also, within the Z2 
slavespin approach one cannot find local moments such 



as an antiferromagnetically ordered phase (AFM), but 
instead obtains a valence bond solid (VBS) phase. In the 
limit A — >• it is obvious that one should find Neel or- 
der instead and that the VBS order is an artifact of the 
specific slave particle approach. 

Regarding the values of U c (e.g. for A = A = 0), one 
should keep in mind that the microscopic U c ~ 4.3 as 
found within QMC^is understimated by slave rotor the- 
ory (U c = 1.68) while it is overestimated by the slave spin 
approach (U c ~ 8) (Fig.[2}. 

IV. STRONG COUPLING LIMIT 

We consider the limit of infinitely strong electron- 
electron interactions. As a result, charge fluctuations 
are frozen out and we obtain a pure spin Hamiltonian at 
half filling. Most importantly, the complex next-nearest 
neighbor spin orbit hoppings result in anisotropic and 
more complicated second neighbor spin exchange terms 
which we analyze in the following. 

A. Kane— Mele spin model 

Taking the limit U — >• 00 of the Kane-Mele-Hubbard 
model Q results in the effective spin model^ 

^km = J\ ^ SiSj + J\ ^ [— Sf Sj — S^Sj + S?Sj] 

(ij) <ij» 

(6) 

where J\ = 4t 2 /U and J\ = 4A 2 /Z7. The second neigh- 
bor exchange term (indicated by <C • ^>) acting merely 
on individual, i.e. triangular sublattices partially frus- 
trates the system. The XF-spin terms prefer ferromag- 
netic order on the individual sublattices which is consis- 
tent with antiferromagnetic order on the original honey- 
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comb lattice; in contrast, the Ising term SfSj favors an- 
tiferromagnetic order on the sublattice competing with 
both the XF-terms and the J\ term. The magnetiza- 
tion, which might point in any direction for J\ = due 
to spin rotational invariance, turns into the XF-plane 
in order to avoid the frustrating part of the J\ t erm^t 
These findings were confirmed within QMC^^^, varia- 
tional cluster approximation (VCA)^, and cluster dy- 
namical mean-field theory (CDMFT) 44 calculations at 
intermediate U/t « 5... 9 and small A. For small J\, 
one can thus employ Hkm to compare other numerical 
approaches against PFFRG method, which we will use 
in the following. 



B. Sodium iridate spin model 

The strong coupling limit of the SIH model is given by 
the spin Hamiltonian 

Hsi = JiJ2 S i S i- J x £ SiSj+2J- x S7S]. (7) 

(ij) <Cr/» 7— links 

Note that the 7-links are the second neighbor links (the 
green, red, and blue lines in Fig. lip). It is struct urally 
similar to the Heisenberg-Kitaev (HK) Hamiltonian^^ 
which has been found to adequately describe the A2lrC>3 
iridates from a spin-orbit Mott picture (A=Na or Li/™. 
Whereas the SI model assumes the nearest-neighbor hy- 
bridization to be trivial and to be essentially given by real 
Ir-Ir hybridization, the kinetic theory underlying the HK 
model more carefully resolves the emergent ter ms fr om a 
multi-orbital Ir-0 cluster superexchange model 4 ^ 4 ^. De- 
pending on the Ir-O-Ir angle, these terms are either more 
or less relevant than the next nearest neighbor exchange 
terms which are considered in the SI model 37 . For the 
links in vertical direction (links with <j z \ one obtains the 
same term as for Hkm, while for the links associated with 
a x one finds +SfS? - Sf S? - Sf SJ and so on. For Hkm 
we have seen that the magnetization turns into the XY- 
plane. Here, however, the term +SfSj — S? 5| — SfS? 
will force the magnetization into the FZ-plane while the 
term -SfS* + S? S? - SfS? favors the XZ-plane and 
so on. Since all the terms (links) are equally distributed 
over the lattice, a priori no plane or direction is preferred. 
As the system sets out to be Neel-ordered for = 0, 
it is conceivable that the competing ordering tendencies 
might at first compensate each other and allow for a per- 
sistent Neel order at small J\. 



V. METHOD 

The PFFRG approacfPH^Eol starts by reformu- 
lating the spin Hamiltonian in terms of a pseudo 
fermion representation of the spin- 1/2 operators S^ = 
V 2 fl a a{3fpi K P =t> Uv> = x, y, z) with fermionic 
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FIG. 3. Characteristic behavior of the flowing (A-dependent) 
susceptibility in a magnetically ordered phase. While the RG 
flow is smooth above some critical value A c ~ 0.45, a numeri- 
cally unstable regime is found below that value. This feature 
signals a magnetic instability which becomes a divergence in 
the thermodynamic limit. The specific case shown here rep- 
resents the largest component of the susceptibility of 'Hsi at 
= 0.4Ji where the system favors antiferromagnetic order. 



operators /f and f± and Pauli-matrices a M . Such a rep- 
resentation enables us to apply Wick's theorem, lead- 
ing to standard Feynman many-body techniques. In this 
pseudofermion language, quantum spin models become 
strongly coupled models with zero fermionic bandwidth 
and finite interaction strength. 

A major advancement of the PFFRGf^is that it allows 
to tackle this situation by providing a systematic scheme 
for the infinite order self-consistent resummations. The 
first conceptual step is the introduction of an infrared fre- 
quency cutoff A in the fermionic propagator. The FRG 
then formulates differential equations for the evolution 
of all m-particle vertex functions under the flow of A^H. 
Hence, one might think of the diagrammatic summations 
as being performed during the RG flow: each discretized 
RG step effectively increases the amount of diagrams in- 
cluded in the approximation. 

To reduce the infinite hierarchy of coupled equations to 
a closed set, a common approach is to restrict oneself to 
one-loop diagrams. The PFFRG extends this approach 
by also including certain two-loop contribution d 35 * 52 1 to 
retain a sufficient backfeeding of self-energy corrections 
to the two-particle vertex evolution. A crucial property 
of the PFFRG is that the the two-particle vertex includes 
both graphs that favor magnetic order and those that fa- 
vor disorder in such a way that the method treats both 
tendencies on equal footing 3 -^. It is the two-particle ver- 
tex which allows to extract magnetic susceptibility as the 
main outcome of the PFFRG. The FRG equations are 
simultaneously solved on the imaginary frequency axis 
and in real space. A numerical solution requires (i) to 
discretize the frequency dependencies and (ii) to limit 
the spatial dependence to a finite cluster, thus keeping 
correlations only up to some maximal length. In our cal- 
culations, the latter typically extends over distances of 
up to 9 lattice spacings corresponding to a correlation 



5 



area (cluster size) of 181 lattice sites of the hexagonal 
lattice. The onset of spontaneous long-range order is sig- 
naled by a sudden breakdown of the smooth RG flow, 
while the existence of a stable solution indicates the ab- 
sence of long-range order. (See Refs . l34l and l35l for further 
technical details.) Fig. [3] shows an example for the char- 
acteristic flow behavior in a magnetically ordered phase. 

VI. RESULTS 

A. Kane— Mele spin model 

From Eq. (J6|, the Kane-Mele spin model reduces to 
an isotropic nearest neighbor spin system in the limit of 
vanishing spin orbit coupling J\ = 0. In this case, the 
system exhibits the standard Neel state on the honey- 
comb lattice. Within our PFFRG approach, this type of 
order is signaled by an instability breakdown in the RG 
flow occurring at the K- and if'-points, i.e. the corners 
of the extended (second) Brillouin zone of the honeycomb 
lattice. (Unless stated otherwise, we plot the susceptibil- 
ity in the second Brillouin zone of the underlying two- 
atomic Bravais lattice because the experimentally con- 
nected unfolded susceptibility has the periodicity of this 
extended zone.) Hence, at an RG scale right before the 
magnetic order sets in, the momentum resolved suscep- 
tibility shows pronounced peaks at the K- and if '-point 
positions. As a consequence of rotational invariance the 
susceptibility profile is identical for all directions of exter- 
nal magnetic fields. Once the spin orbit interaction J\ is 
switched on, the situation changes considerably as shown 
in Fig. [4] While the susceptibility peaks for an external 
field in ^-direction (or ^/-direction) become even sharper 
as compared to J\ = 0, the peaks in the z-component 
drop drastically. Already at small J\ = 0.1 this effect is 
rather pronounced, which evidences that for finite Ja, the 
spins favor the x-y plane. (We set J\ = 1 in this section.) 
With increasing Ja, more weight of the z-susceptibility is 
transferred to the x- and ^/-components of x- For strong 
enough Ja, the remnant magnetic fluctuations in \ z are 
not of antiferromagnetic type anymore, which can be seen 
in Fig. [4] for J\ = 0.5 showing small maxima at M-point 
positions. We do not observe any particular phase transi- 
tion at A > 0. In particular, the magnetic order persists 
in the whole parameter space. The frustration gener- 
ated by the J\S* 5| -terms has little effect because the 
spins can circumvent this frustration by avoiding the z- 
axes. With increasing Ja, the two sublattices become 
effectively decoupled such that in the limit Ja — >• 00 both 
sublattices exhibit xy ferromagnetic order independently. 

B. Sodium iridate spin model 

As for the KM spin model in the previous section, the 
SI spin model becomes a simple isotropic nearest neigh- 
bor spin system in the limit = and hence shows 




FIG. 4. Magnetic susceptibilities at the critical scale A = A c 
for various values of Ja (Ji = 1) in the Kane-Mele spin 
model, resolved for in plane (x,y) and out of plane (z). Top 
row: x x (k) (left panel) and x z {k) (right panel) for Ja = 0.1. 
Bottom row: x x (M) = X V (k) (left panel) and x*(&) (right 
panel) for J\ — 0.5. The susceptibility weight along z signif- 
icantly decreases for large Ja. For higher Ja, the remainder 
z-susceptibility deviates from the Neel AFM structure. 



Neel order (upper left plot in Fig. [5]). For finite but not 
too large Js^, the antiferromagnetic order persists, i.e. 
the position of the ordering peaks in the susceptibility 
remains unchanged (J^ = 0.5 in Fig. |5|. As the suscep- 
tibility looses its sixfold rotation symmetry for finite J^, 
this manifests in the deformation of the ordering peaks as 
compared to J^ =0. Note that due to the special connec- 
tion between lattice directions and spin directions in the 
SI spin model, the a>, y- and ^-components of the suscep- 
tibility transform into each other under /c-space rotations 
of 120° in clockwise direction. Fig. [5] illustrates x z which 
preserves the symmetries k x — >• — k x and k y — >• —k y . Note 
that regardless of the particular phase, the value of the 
susceptibility at the six K- and if '-points must always be 
equal. This results from the fact that the three if-points 
(or if'-points) are related by reciprocal lattice vectors 
among each other. Furthermore, since the two sublat- 
tices are equivalent, the K- and if '-points are likewise 
degenerate. 

An interesting observation can be made regarding the 
orientation of the antiferromagnetic order. Due to the 
equivalence of the x- : y- and z-direction in spin space, 
the magnetic order can point in each of these directions 
without any preference. Even though SU(2) symmetry is 
explicitly broken, the rotational symmetry of the suscep- 
tibility prevails: consider a magnetic field B = vJ? point- 
ing in some direction v = 3, y z v v e v with |v| = 1. 
The corresponding susceptibility x v ' •> i- e - the linear re- 
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FIG. 5. Magnetic susceptibilities of the SI spin model for various values of (Ji = 1). All susceptibilities shown refer to 
a magnetic field in z-direction. The x- and ^-components of th e susc eptibility are obtained by /c-space rotations of 120° in 
clockwise or counterclockwise direction, respectively. (See Section VI B for more details.) While the Neel peaks initially persist 
for finite J^, the peaks start to move due to the onset of incommensurability (Fig. |6|. For large J^, new suceptibility peaks 
emerge which link to the change of unit cell structure of magnetic order. 



sponse to such a perturbation is defined as 



X 



dB \b^o dB 

J2 9 ^=x,y,z dBj, 
H'=x,y,z 



dB 



(8) 



= E v ^'v 



where 



fi,/i =x,y,z 



dM u 



and M is the magnetization. 



dB^ IB^O 

Since x M/i cannot develop any off-diagonal elements be- 
fore reaching the magnetic instability in the RG flowj^ 
we have x MM = Sllll'X^- It follows that 



(j,=x,y,z 

Since x x — X y — X z a ^ all ) -points, we obtain 



E 



(9) 



(10) 



Hence, in linear response the low energy physics of the 
system is rotationally symmetric and the antiferromag- 
netic order can point in any direction. This is a con- 
sequence of the Neel order residing at high-symmetry 
points of the Brillouin zone as well as the special con- 
nection between lattice directions and spin directions in 
the SI spin model. However, this argument does not hold 
for spin fluctuations away from the K- or if '-points. For 
fluctuations at arbitrary momentum, a certain direction 
will generally be preferred. 



>, 5 2n 




FIG. 6. Dependence of the ordering vector Q y on in Hsi- 
The inset illustrates the evolution of the ordering peaks in the 
Brillouin zone (thick hexagon: second Brillouin zone, thin 
hexagon: first Brillouin zone; see also Fig. |5|. In the limit 
— » oo, the system converges again towards a commensurate 
ordering vector. 



As increases, the deformation of the ordering peaks 
at the K- and if '-points becomes more pronounced. At 
some coupling w 0.53, the peaks split and the new 
maxima move along the k y -direction (Fig. [6]). These peak 
positions indicate a phase transition to a spiral phase 
with incommensurate order. It is important to note, how- 
ever, that magnetic order persists in the whole parameter 
regime around the transition and we find no magnetically 
disordered phase. This can be seen from the behavior of 
the RG flow which always exhibits a characteristic insta- 
bility breakdown. To demonstrate the evolution of the 
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ordering vector in the spiral phase, Fig. [6] shows the peak 
position as function of J^. Note that the k x -component 
of the peak position is constant in J^. With increas- 
ing J^, the peaks move continuously towards the points 
Qoo = ±§ ^7=) which lie at two third of the distance 

between the ) -points and the k x -ax.is (Fig. [H|). Again, 
with increasing there is no sign of any non-magnetic 
phase. 

The system at infinite spin-orbit coupling is of particu- 
lar interest, as this case represents a model with Kitaev- 
like interactions on the triangular lattice. As goes to 
infinity, the system is effectively described by decoupled 
triangular sublattices. Hence, already the first Brillouin 
zone, i.e. the Brillouin zone of a triangular sublattice, 
contains the full information about x m &-space. The 
susceptibility then becomes periodic with respect to this 
smaller zone. Such a change of periodicity can be seen in 
Fig. [5] at large where new peaks at k x = emerge. In 
the limit — >• 00 these new peaks reach the same height 
as the ones at and finally become identical to them, 
indicating the new periodicity in /c-space. Fig. [7^i shows 
the susceptibility in the first Brillouin zone of the trian- 
gular sublattice in this limit. From the peak positions, 
one can easily derive the corresponding spin pattern in 
real space. On each triangular sublattice the wave vec- 
tor is half the one of the 120°-Neel order residing at the 
corners of the first Brillouin zone. The unit cell contains 
6x6 lattice sites as compared to the 3x3 unit cell of the 
120°-Neel order. Hence, the order is commensurate and 
the local magnetic moments along a lattice direction are 
modulated with a periodicity of 6 sites. Taking into ac- 
count both sublattices of the honeycomb lattice, we end 
up with a unit cell containing 72 sites. 

Our numerical conclusions for the SI spin model in the 
limit — >• 00 can also be reconciled with an analytical 
argument. Performing a transformation in spin space, 
Si — >• the system at this point can be mapped to 
an SU(2) invariant antiferromagnetic Heisenberg model 
on the triangular lattice, Hsi = Y^ij Si^j- For this map- 
ping, we divide the triangular lattice into four sublattices 
denoted by "xy" , "xz" and "yz" , each with a dou- 
bled lattice constant (Fig. JtJd). The relation between S$ 
and Si depends on the sublattice, 



i G " • " : 


Q _ ( QX QV qz\ 
°i — \°i i°i 5°i ) 1 




i e "xy" : 


Si — ( — Sf , — Sf , S^) , 




i £ "xz v : 


S$ = (—Sf, Sf, —Sf) , 




i e "yz" : 


Q f QX qV qz\ 


(ii) 



e.i., while on sublattice "•" the spins remain unchanged, 
on the sublattice "xy" the x- and ^/-components of the 
spin operator acquire a minus sign, and so on (a similar 
mapping for the Heisenberg-Kitaev model at a = 0.5 
is described in Ref. 08]). Since the antiferromagnetic 
Heisenberg model on the triangular lattice exhibits mag- 
netic order via the 120°-Neel stated it follows that the 
SI spin model at — » 00 is likewise magnetically or- 




(b) 




FIG. 7. The SI spin model at — >> 00: (a) Magnetic suscep- 
tibility displayed in the first Brillouin zone of the triangular 
sublattice. The two ordering peaks correspond to the peaks 
in Fig. [5] which emerge at > 5 and k x — 0. In the limit 

— > 00, these maxima reach the same hight as the ones at 
Qoo = (±^f,±§3|). (b) Mapping of the SI spin model at 

— >• 00 to the antiferromagnetic Heisenberg model on the 
triangular lattice: The lattice is divided into four sublattices 
denoted by "xy", "xz" and "yz". As shown in Eq. (11) 
the transformation from Si to Si depends on the sublattice 
where i resides. The exchange couplings follow the convention 
shown in Fig. [l] 



dered. The corresponding spin pattern in real space can 
be found by applying the inverse of the above spin trans- 
formation to the 120°-Neel state: Since the structure of 
the spin rotations (Fig. [7)3) has a periodicity of two lat- 
tice sites in each lattice direction, the 3x3 unit cell of 
the 120°-Neel order transforms back into a 6 x 6 unit cell, 
as found within our PFFRG calculations. 



VII. DISCUSSION 

In view of our results for the SI model, we specu- 
late about the implications for the phase diagram at in- 
termediate U. We find the incommensurate phase for 
J-JJx > 0.53 where J~ x = 4X 2 /U and J x = 4t 2 //7, imply- 
ing a transition at j w 0.73 for large U. In Fig. [ sj w e 
have replotted the phase diagram of Riiegg and FietlP^ in 
a slightly modified way. Our reasoning is the following: 
since the spin model corresponds to U — >• 00, we extrap- 
olated the phase boundary between "VBS (AFM)" and 
"QSH*" of the phase diagram in Ref. [33] to larger U. As 
the phase transition occurs for sufficiently large U ap- 
proximately at j ~ 0.73 (Fig. 8), we speculate that the 
observed phase transition from Neel to spiral order in the 
spin model is a remnant of the phase transition into the 
QSH* phase at intermediate U . Within this scenario, 
the QSH* phase would transform into spiral magnetic 
order in the limit U —> 00. We note, however, that this 
necessarily implies that with increasing U the gap of the 
QSH* phase closes at some point to form the Goldstone 
mode of the spiral order. In principle, the gap closure 
can occur at finite U (which would imply an additional 
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FIG. 8. (color online) Schematic phase diagram for the SI 
Hubbard model as conjectured from our strong coupling re- 
sults of a Neel to spiral transition at = 0.53 (upper x 
scale), corresponding to X/t ~ 0.73 (lower x scale) . While 
we cannot ultimately assess the nature of a possible interme- 
diate phase aside from Tl and AFM for strong interaction 
and strong spin-orbit coupling (yellow phase) , its existence is 
likely. 



doubled semion spin liquid^. The alternative scenario — 
assuming that a QSH*-type phase does not exist — would 
still require an additional phase comparable to the yel- 
low phase of the schematic phase diagram in Fig.|8|; in 
this case, the additional phase would most likely be a 
magnetically ordered phase {e.g. spiral order). Whether 
or not this phase is a topological liquid or just another 
magnetically ordered phase, we conjecture that in either 
case an additional phase of some kind should be present. 

In summary, we find that the physics of the SI spin 
model is much richer as compared to the KM spin model. 
This can be traced back to the different spin symme- 
tries in both systems. The broken axial symmetry in 
the SI spin model prevents the spins from forming pla- 
nar antiferromagnetic order and eventually leads to the 
emergence of a spiral phase. Note that this phase does 
not have any analogue in similar models such as the 
Heisenberg-Kitaev model As such, we have identified 
multi-directional spin orbit terms to be an interesting 
way to create new spin phases in the infinite U limit and 
possibly even more exotic phases at intermediate U when 
charge fluctuations enter the picture. 

To give another direction of further investigation, it 
will be interesting to study the KM model in the presence 
of Rashba spin orbit coupling 



(ij) a(3 



phase boundary in Fig. [8]) or at U —> oo. In the latter 
scenario, the QSH* phase could extend up to U —> oo. 
Furthermore, one should keep in mind that the QSH* 
phase as found in Ref. [33] might not be the only topo- 
logical liquid candidate with similar properties, such as a 



The Rashba term breaks the remaining U(l) symmetry of 
the electron spin to Z2, and also affects the z — >• — z mir- 
ror symmetry as well as particle-hole symmetry. Taking 
into account the Rashba spin orbit coupling results in a 
more complicated spin Hamiltonian with some terms be- 
ing of Dzyaloshinskii-Moriya type. The full Hamiltonian 
is given by 



^KMR = J\ ^ [—SfSj — S^Sj + S?Sj] 

[(Ji + Wf S* + (Ji - Jr)(S?S] + S?S?) - y/jJniSfS* - SfS?) 



S 1 — links 

■ E 

5^— links 

- E 

£3 — links 



b i/3 ^ 



QX 



1 QX 



"2 w 



^Sfr, and S* -^Sf. - is* 



S" : and S* 



2 i H 



if 3 



i qy 

2°i/3 



(12) 



where the third line in (|12| is obtained from the second 
one by replacing S^- by —1/2S^- — y/^/2S^. and so 
on. The different links denoted by Si (i = 1, 2, 3) are the 
three nearest neighor vectors of the honeycomb lattice. 
We expect the J\-Jr phase diagram to be interesting 
and to host some additional phases, which we defer to a 



future publication. 
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VIII. CONCLUSION 



We have investigated the strong coupling limit of Hub- 
bard models of topological honeycomb band structures. 
We have considered two band structures both classified 
as two-dimensional Z2 topological insulators, where only 
the Kane-Mele spin model as opposed to the sodium iri- 
date model preserves axial spin symmetry. For the for- 
mer model at infinite coupling, the magnetism tends to 
form XY antiferromagnetic order already at very small 
spin orbit couplings. This way the spins manage to avoid 
the frustration induced by the spin-orbit anisotropic spin 
terms. As a consequence, frustration effectively plays no 
role in the KM spin model. The physical scenario is very 
different for the sodium iridate model with generalized 
spin orbit couplings and hence broken axial spin symme- 
try. There, the spins cannot form XY , XZ or YZ order. 
As a result, the magnetic phase formation in the strong 
coupling limit exhibits a commensurate to incommensu- 
rate Neel to spiral transition at ~ 0.53. In the limit 
of infinite spin orbit coupling, the model converges to a 
commensurate magnetic state with a 6 x 6-site unit cell 
on each of the two decoupled triangular sublattices of the 
underlying honeycomb model. The emergence of the spi- 
ral phase in the infinite U limit leads us to conjecture that 



aside from the topological band insulator regime and the 
antiferromagnetic phase, a third phase should exist at fi- 
nite U and finite spin orbit coupling. In this respect, our 
results are not inconsistent with the existence of a frac- 
tionalized QSH* phase as proposed in Ref.[33j Whatever 
this phase will eventually turn out to be, we find that 
the breaking of axial spin symmetry is generally vital to 
the emergence of new phases and an enriched diversity 
of magnetic phases in interacting topological honeycomb 
band structures. 

Note added. Recently, the classical Kitaev-Heisenberg 
model has been studied on the triangular lattice within 
Monte Carlo approaches. Interestingly, spiral-like order 
has been found to be associated with a Z 2 vortex lat- 
ticed 
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